Fast simulation of stabilizer circuits using a graph state representation 



Simon Anders 1 'Q and Hans J. Briegel 1 ' 2 

1 Institut fiir Theoretische Physik, Universitdt Innsbruck, Innsbruck, Austria 
2 Institut fiir Quantenoptik und Quanteninformation der Osterreichischen Akademie der Wissenschaften, Innsbruck, Austria 

(Dated: December, 2005 (v2)) 

According to the Gottesman-Knill theorem, a class of quantum circuits, namely the so-called 
stabilizer circuits, can be simulated efficiently on a classical computer. We introduce a new algorithm 
for this task, which is based on the graph-state formalism. It shows significant improvement in 
comparison to an existing algorithm, given by Gottesman and Aaronson, in terms of speed and of 
the number of qubits the simulator can handle. We also present an implementation. 
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I. INTRODUCTION 

Protocols in quantum information science often use en- 
tangled states of a large number of qubits. A major chal- 
lenge in the development of such protocols is to actually 
test them using a classical computer. This is because 
a straight-forward simulation is typically exponentially 
slow and hence intractable. Fortunately, the Gottesman- 
Knill theorem (0,0) states that an important subclass 
of quantum circuits can be simulated efficiently, namely 
so-called stabilizer circuits. These are circuits that use 
only gates from a restricted subset, the so-called Clif- 
ford group. Many techniques in quantum information 
use only Clifford gates, most importantly the standard 
algorithms for entanglement purification EL 1 1 , 0, IE 
and for quantum error correction |E llOl lll| . Hence, 
if one wishes to study such networks, one can simulate 
them numerically. 

The usual proof of the Gottesman-Knill theorem (as 
stated e. g. in 0) contains an algorithm that can carry 
out this task in time 0(N 3 ), where N is the number of 
qubits. Especially for the applications just mentioned, 
one is interested in a large N: For entanglement pu- 
rification one might want to study large ensembles of 
states, and for quantum error correction concatenations 
of codes. The cubic scaling renders this extremely time- 
consuming, and a more efficient algorithm should be of 
great use. 

Recently, Aaronson and Gottesman presented such an 
algorithm (and an implementation of it) in Ref. |12| . 
whose time and space requirements scale only quadrati- 
cally with the number of qubits. In the present paper, we 
further improve on this by presenting an algorithm that 
for typical applications only requires time and space of 
0(N log N). While Aaronson and Gottesman's simula- 
tor, when used on an ordinary desktop computer, can 
simulate already systems of several thousands of qubits 
in a reasonable time, we have used our simulator for over 
a million of qubits. This provides a valuable tool for in- 
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vestigating complex protocols such as our study of multi- 
party entanglement purification protocols in Ref. [l3| . 

The crucial new ingredient is the use of so-called graph 
states. Graph states have been introduced in [l4| for the 
study of entanglement properties of certain multi-qubit 
systems; they were used as starting point for the one- 
way quantum computer (i. e., measurement-base quan- 
tum computing) |lq |. and found to be suited to give 
a graphical description of CSS codes (for quantum er- 
ror correction) Graph states take their name from 
the concept of graphs in mathematics: Each qubit cor- 
responds to a vertex of the graph, and the graph's edges 
indicate which qubits have interacted (see below for de- 
tails). 

There is an intimate correspondence between stabilizer 
states (the class of states that can appear in a stabilizer 
circuit) and graph states: Not only is every graph state 
a stabilizer state, but also every stabilizer state is equiv- 
alent to a graph state in the following sense: Any sta- 
bilizer state can be transformed to a graph state by ap- 
plying a tensor product of local Clifford (LC) operations 
[171 lia ITgj . We shall call these local Clifford operators 
the vertex operators (VOPs). 

To represent a stabilizer state in computer memory, 
one stores its tableau of stabilizer operators, which is 
an N x N matrix of Pauli operators and hence takes 
space of order 0(N 2 ) (see below for details). Gottesman 
and Aaronson's simulator extends this matrix by another 
matrix of the same size (which they call the destabilizer 
tableau), so that their simulator has space complexity 
0(N 2 ). A graph state, on the other hand, is described 
by a mathematical graph, which, for reasons argued later, 
only needs space of O(NlogN) in typical applications. 
Hence, much larger systems can be represented in mem- 
ory, if one describes them as graph states, supplemented 
with the list of VOPs. However, we also need efficient 
ways to calculate how this representation changes, when 
the represented state is measured or undergoes a Clifford 
gate application. The effect of measurements has been 
extensively studied in [2E|, and gate application is what 
we will study in this paper, so that we can then assemble 
both to a simulation algorithm. 

This paper is organized as follows: We first review the 
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stabilizer formalism, the Gottesman-Knill theorem, and 
the graph state formalism in Section |H| There, we will 
also explain our representation in detail. Section IlIII cx- 
plains how the state representation changes when Clif- 
ford gates are applied. This is the main result and the 
most technical part of the paper. For the simulation of 
measurements, we can rely on the studies of Ref. |20j . 
which are reviewed and applied for our purpose in Sec- 
tion El Having exposed all parts of the simulator algo- 
rithm, we continue by presenting our implementation of 
it. A reader who only wishes to use our simulator and 
is not interested in its internals may want to read only 
this section. Section IVII assesses the time requirements 
of the algorithm's components described in Sections IIIII 
and IIVI in order to prove our claim of superior scaling of 
performance. We finish with a conclusion fSection |VII|) . 



II. STABILIZER AND GRAPH STATES 

We start by explaining the concepts mentioned in the 
introduction in a formal manner. 

Definition 1. The Clifford group Cjy on N qubits is de- 
fined as the normalizer of the Pauli group Vjq : 

C N = {U e SU(2 N )\UPU^ eV N VP e V N ] , 

V N = {±l,±i}-{I,X,Y,Z}® N , (1) 

where I is the identity and X , Y , and Z are the usual 
Pauli matrices. 

The Clifford group can be generated by three elemen- 
tary gates (see e. g. 0): the Hadamard gate H, the f 
phase rotation S, and a two-qubit gate, either the con- 
trolled NOT gate AX, or the controlled phase gate AZ: 
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The significance of the Clifford group is due to the 
Gottesman-Knill theorem (0, see also 0): 

Theorem 1. A quantum circuit using only the follow- 
ing elements (called a stabilizer circuit) can be simulated 
efficiently on a classical computer: 

• preparation of qubits in computational basis states 

• quantum gates from the Clifford group 

• measurements in the computational basis 



Definition 2. An N-qubit state is called a stabilizer 
state if it is the unique eigenstate with eigenvalue +1 of 
N commuting multi-local Pauli operators P a ( called the 
stabilizer generators) : 

P*\tl)) = W), P a eV N , a=l,...,N 

(These N operators generate an Abelian group, the 
stabilizer, of 2 N Pauli operators that all satisfy this sta- 
bilization equation.) 

Computational basis states are stabilizer states. Fur- 
thermore, if a Clifford gate U acts on a stabilizer state 
the new state U\ip) is a stabilizer state with genera- 
tors UPiW G Vn- Hence, the state in a stabilizer circuit 
can always be described by the stabilizer tableau, which 
is a matrix of N x N operators from {I, X, Y, Z} (where 
each row is preceded by a sign factor). The effect of an 
n-qubit gate can then be determined by updating nN 
elements of the matrix, which is an efficient procedure. 

Instead of on the stabilizer tableau, we shall base our 
state representation on graph states: 

Definition 3. An N-qubit graph state \G) is a quantum 
state associated with a mathematical graph G — (V,E), 
whose \ V\ — N vertices correspond to the N qubits, while 
the edges E describe quantum correlations, in the sense 
that \G) is the unique state satisfying the N eigenvalue 
equations 

K { « ] \G) = \G), aeV, 

withK^=a^ 11 a^=:X a J] Z b , (3) 

fc£ngbha bGngbha 

where ngbha := \b I \a, b} € E} is the set of vertices ad- 
jacent to a flA . \lA. 

The following theorem states that the edges of the 
graph can be associated with phase gate interactions be- 
tween the corresponding qubits: 

Theorem 2. // one starts with the state \+)® N — 
Yia^v H a \00 ■■■ 0) one can easily construct \G) by apply- 
ing AZ on all pairs of neighboring qubits: 



\G) = 



[I AZ ab 
K {a,b}eE 



n ^ i°) c 



(4) 



\a£V 



(Proof: Insert Eq. © into Eq. © HJ.) 

As the operators Kq belong to the Pauli group, all 
graph states are stabilizer states, and so are the states 
which we get by applying local Clifford operators C G C\ 
to \G) . For such states, we introduce the notation 



\G;C_) :— |G; Ci, C2, . . . , Cm) 



iQIG) (5) 



The proof of the theorem is simple after one introduces 
the notion of stabilizer states [2l|: 



It has been shown that all stabilizer states can be 
brought into this form [TtI ITsl Il9|. i. e. any stabilizer 
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state is LC-equivalent to a graph state. (We call two 
states LC-equivalent if one can be transformed into the 
other by applying a tensor product of local Clifford op- 
erators.) Finding the graph state that is LC-equivalent 
to a stabilizer state given by a tableau can be done by a 
sort of Gaussian elimination as explained in |17| . 

This is what we shall use to represent the current quan- 
tum state in the memory of our simulator. Fig. ^ shows 
for an example state the tableau representation that is 
usually employed (and also used by CHP, albeit in a mod- 
ified form) and our representation. The tableau repre- 
sentation requires space of order 0(N 2 ). We store the 
graph in adjacency list form (i. e., for each vertex, a list 
of its neighbors is stored), which needs space of order 
0{Nd), where d is the average vertex degree (number 
of neighbors) in the graph. We also store a list of the 
N local Clifford operators C\ , . . . , Cn , which transform 
the graph state \G) into the stabilizer state \G;C). We 
call these operators the vertex operators (VOPs). As 
there are only 24 elements in the local Clifford group, 
each VOP is represented as a number in 0, . . . , 23. The 
scheme to enumerate the 24 operators will be described 
in |22j. Note that we can disregard global phases of the 
VOPs as they only lead to a global phase of the full state 
of the simulator. 

As we shall see later, we may typically assume that d = 
0{\ogN). Hence, our representation needs considerably 
less space in memory than a tableau, namely 0(N log N), 
including O(N) for the VOP list. 

The Gaussian elimination needed to transform a sta- 
bilizer tableau to its graph state representation is slow 
(time complexity 0(N 3 )), and so we should better not 
use it in our simulator. But usually, one starts with the 
initial state jO)®^, and if we write this state already in 
graph state form, the tableau representation is never used 
at all. 

From Eq. it is clear that the initial state can be 
written as a graph with no edges and Hadamard gates 
acting on all vertices: 

\0)® N = |({1,..., N}, {});H,...,H). 



III. GATES 

When the simulator is asked to simulate a Clifford gate, 
the current stabilizer state is changed and its graph repre- 
sentation has to be updated to correctly reflect the action 
of the gate. How to do this, is the main technical result 
of this paper. 
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FIG. 1: A stabilizer state \ip) represented in different ways: 

(a) as stabilizer tableau, i. e. the state is stabilized by the 
group of Pauli operators generated by the operators in the 4 
rows. This representation needs space 0(N 2 ) for N qubits. 

(b) , (c) as LC-equivalence to a graph state, (b) shows the 
graph, with the VOPs given by their decomposition into the 
group generators {H, S}. (c) is the data structure that repre- 
sents (b) in our algorithm. The VOPs are now specified using 
numbers between and 23 (which enumerate the \Ci\ = 24 
LC operators). Here, we need space O(Nd), where d is the 
average vertex degree, i. e. the average length of the adja- 
cency lists. Writing G for the graph in (b), we can use the 
notation of Eq. Q and write \%j>) = \G;H, I, HS, S). 



2. Two-qubit gates 

It is sufficient if the simulator is capable to simulate 
a single multi-qubit gate: As the entire Clifford group 
is generated, e. g., by H, S, and AZ, all gates can be 
constructed by concatenating these. We chose to imple- 
ment AZ, the phase gate, as this is (because of its role 
in Eq. (0}) most natural for the graph-state formalism. 

In the following discussion, the two qubits onto which 
the phase gate acts, are called the operand vertices and 
denoted with a and b. All other qubits are called non- 
operand vertices and denoted c, d, . . . . 

To solve the task, we have to distinguish several cases. 

Case 1. The VOPs of both operand vertices are in 
Z, where Z := {I, Z, S, S^} denotes the set of those four 
local Clifford operators that commute with AZ (the other 
20 operators do not). In this case, applying the phase 
gate is simple: We use the fact that (due to Eq. (QJ) 
applying a phase gate on a graph state just toggles an 
edge: 

AZ ab \(V,E)) = \(V,E A {{a,b}})), 



1. Single- qubit gates 

In the graph representation, applying local (single- 
qubit) Clifford gates becomes trivial: if C £ C\ is applied 
to qubit a, we replace this qubit's VOP C a by CC a - 



where A denotes the symmetric set difference 4a6:= 
(Ai)B)\(Ar\B), i. e. the edge {a, b} is added to the graph 
if is was not present before, otherwise it is removed. 

Case 2. The VOP of at least one of the operand ver- 
tices is not in Z. In this case, just toggling the edge 
is not allowed because the AZ a t cannot be moved past 
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the non-2 VOP. But there is a way to change the VOPs 
without changing the state, which works in the following 



'iX to C a . This factor 



Both operand vertices have non- 
Here, the following operation will 



Sub-case 2.2. 

operand neighbors. 
help: 

Definition 4. The operation of local complementation 
about a vertex a of a graph G = (V,E), denoted L a , 
is the operation that inverts the subgraph induced by the 
neighborhood of v: 

L a (V, E) = (V, E A {{b, c}\b, c e ngbha}) 



This operation transforms the state into a local- 
Clifford equivalent one, as the following theorem, taken 
from [I3,[2(|, asserts: 

Theorem 3. Applying the local complementation L a 
onto a graph G yields a state \L a G) = U\G), with the 
multi-local unitary 



U = y/-iX a ]"[ 



bfEngbh a 

Note that the operator yfiZ is related to the phase op- 
erator S of Eq. ©: ViZ = e l T S\ and sfiX = <J=iX ] = 

An obvious consequence of Theorem |3| is the following. 
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Corollary 1. A state \G;C_) is invariant under applica- 
tion of L a to G, followed by an updating of C according 
to 



CbViX for b = a 
Cb^ { CbV-iZ for b E ngbha 
Cb otherwise 



(6) 



Now note that the local Clifford group is generated not 

>iZ, 



only by S and H but also by \/—iX and ViZ, the Her- 
mitian adjoints of the operators right-multiplied to the 
VOPs in Eq. © . Our simulator has a look-up table that 
spells out every local Clifford operator as a product of 
-as it turns out, at most 5- of these two operators, times 
a disregarded global phase. For example, the table's line 
for H reads: 



H oc \Z—iXViZViZ\/iZV—iX. 



(7) 



This allows us now to reduce the VOP C a of any non- 
isolated vertex a to the identity I by proceeding as fol- 
lows: The decomposition of C a taken from the look-up 
table is read from right to left. When a factor \J —iX is 
read we do a local complementation about a. This does 
not change the state if the correction of Eq. © is applied, 



which right-multiplies a factor 
\fiX cancels with the factor 



-iX at the right-hand end 
of C a 's decomposition, so that we now have a VOP with 
a shorter decomposition. 

If the right-most operator of the decomposition is \fiZ 
we do a local complementation about an arbitrarily cho- 
sen neighbor of a, called a's "swapping partner" . Now, 
the correction operation will lead to a factor S being 
right-multiplied to C a , again shortening the decomposi- 
tion. 

Note that a local complementation about a never 
changes the edges incident on a and hence, if a was non- 
isolated in the beginning of the procedure, it will stay 
so. This is important, as only a non-isolated vertex can 
have a swapping partner. Hence, the procedure can be 
iterated, and (as the decompositions have a maximum 
length of 5) after at most 5 iterations, we are left with 
the identity I as VOP. 

We apply the described "VOP reduction procedure" to 
both operand vertices. After that, both vertices are the 
identity, and we can proceed as in Case 1. 

One might wonder, however, whether the use of the 
VOP reduction procedure on the second operand vertex 
b spoils the reduction of the VOP of the first operand a. 
After all, a could be a neighbor of b or of the swapping 
partner c of b. Then, if a local complementation Lb or 
L c is performed, the compensation according to Eq. © 
changes the neighborhood of b and c (which include a). 
But note that a neighbor of the inversion center only gets 
a factor V 1 —iZ oc S'. As generates Z, this means that 
after the reduction of b, the VOP of a might be no longer 
the identity but it is still an element of Z, and we are 
allowed to go on with Case 1. 

But what happens, if one of the vertices does not have 
a non-operand neighbor, that could serve as swapping 
partner? This is the next Sub-case. 

Sub-case 2.2. At least one of the operand vertices is 
isolated or only connected to the other operand vertex. 
We first assume that the other vertex is non-connected 
in the same sense: 

Sub-sub-case 2.2.1. Both operand vertices are either 
completely isolated, or only connected with each other. 
Then, we can ignore all other vertices and have to study 
only a finite, rather small number of possible states. 

Let us denote by • • the 2-vertex graph with no edges, 
and by •— • the 2-vertex graph with one edge. There are 
only very few possible 2-qubit stabilizer states, namely 
those in 



S 2 :={\G;C 1 ,C 2 ) |Ge{ 



},Ci,C 2 ed}. (8) 



Of course, many of the assignments in the r.h.s describe 
the same state, such that < 2 • 24 2 . Remember that 
the phase gate AZi$ (being a Clifford operator) maps S 2 
bijectively onto itself. 

The function table of AZ 1>2 \s 2 ■ \G;C 1 ,C 2 ) h-> 
\G';C' 1: C2) can easily be computed in advance (we did 
it with Mathematica) and hard-coded into the simulator 
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as a look-up table. This table contains 2 • 24 2 lines such 
as 

|« •,C[i3],C[ 2 ]) H-> |»-»,C [0 ],C [2 ]), (9) 

where the Cta (i = 0, . . . , 23) are the Clifford operators in 
the enumeration detailed in 0| (e. g. C[ ] = ^, CW = V). 

Note that many of the assignments to C\ and C 2 in 
Eq. JSJ describe the same state. Hence, we have a choice 
in the operators C[, C' 2 with which we represent the re- 
sults of the phase gate in the look-up table. It turns out 
(by inspection of all the possibilities) that we can always 
choose the operators such that the following constraint 
is fulfilled: 

Constraint 1. If C\(C2) € Z, choose C' 1 ,C' 2 such 
that again C[(C 2 ) E Z. 

The use of this will become clear soon. 

Sub-case 2.2.2. We are left with one last case, 
namely that one vertex, let it be a, is connected with non- 
operand neighbors, but the other vertex b is not, i. e. has 
either no neighbors or only a as neighbor. Then, we pro- 
ceed as follows: We use iterated local complementations 
to reduce C a to /. After that, we may use the look-up 
table as in Sub-sub-case 2.2.1. That this is allowed even 
though a is connected to a non-operand vertex is shown 
in the following: First note that the state after the re- 
duction of C a to / can be written (following Eq. (JSJ) as 

\{V,E)-C) = l[C c J] AZ cd \+ + ---+) 

c£V {c,d}GE 

= n ° c n Az cd c b (Az Qb ) c i++...+> 

ce {c,d}e 

V\{aM E\{{a,b}} v , 

^ Cb Zd Az ab ' =\+)® N ~ 2 ®M ab 

commute with this with |i^)e<S 2 

(*) 

(where £ = 0,1 indicates whether {a, b} E E). Observe 
that Cb has been moved past the operators AZ c d- This 
is allowed because none of the AZ ca - acts on b 

We now apply AZ ao to this state. AZ ao can be moved 
through all the phase gates and vertex operators above 
the left brace so that it stands right in front of the 1S2 
state \f) ab which is separated from the rest. Thus, the 
table @ from Sub-sub-case 2.2.1 may be used. (This 
would not be the case if, in the state above the brace 
marked with "(*)", the two operand vertices were still 
entangled with other qubits.) The table look-up will give 
new operators C' a , C' b and a new (' , so that the new state 
has the following form: 

AZ ab \(V,E);C) = 

n c c n Az cd c a c b (Az ab f !++■■■+) 

ce {cd}e 

V\{a.b} E\{{a.b}} 

(10) 



For this to be a state in our usual \G;C) form J5j|, 
the two operators C' a and C' b have to moved to the left, 
through the AZ c d- For C b , this is no problem, as b was 
assumed to be either isolated or connected only to a, 
so that C b commutes with Yi{ c d}eE\{{a b}} ^cd, as the 
latter operator does not act on b. The vertex a, however, 
has connections to non-operand neighbors, so that some 
of the AZ c d act on it. We may move it only if C' a E Z (as 
this means that it commutes with AZ). Luckily, due to 
Constraint 1 imposed above, we can be sure that C' a E Z, 
because C a = I S Z. 

Listing 1 shows in pseudo-code how these results can 
be used to actually implement the controlled phase gate 
AZ. 



IV. MEASUREMENTS 

In a stabilizer circuit, the simulator may be asked at 
any point to simulate the measurement of a qubit in the 
computational basis. How the outcome of the measure- 
ment is determined, and how the graph representation 
has to be updated in order to then represent the post- 
measurement state will be explained in the following. 

To measure a qubit a of a state \G,C_) in the computa- 
tional basis means to measure the qubit in the underlying 
graph state \G) in one of the 3 Pauli bases. Writing the 
measurement outcome as £, this means: 

i±£2f*| G ,fi) = ( n c t )l±i^aAG) 

\bev\{a} J 

J n a W + <-'f^ | G) (u, 

\6enW / 

As C a is a Clifford operator, P a := C\Z a C a E 
{X a , Y a , Z a , —X a , -Y a , -Z a }. Thus, in order to measure 
qubit a of \G,C) in the computational basis, we measure 
the observable P a on \G). Note that in case that P a is 
the negative of a Pauli operator, the measurement result 
£ to be reported by the simulator is the complement of 

C, the result given by the X , Y or Z measurement on the 
underlying graph state \G). 

How is the graph G changed and how do the vertex op- 
erators have to be modified if the measurement I± 2 Pa \ G) 
is carried out? This has been worked out in detail in 
Ref. po| . which we now briefly review for the present 
purpose. 

The simplest case is that of P = ±Z. Here, the state 
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cphase (vertex a, vertex b): 
if ngbha\{6} / {}: 

remove_V0P (a, 6) 
end if 

if ngbhb\{a} {}: 

remove_V0P (6, a) 
end if 

[It may happen that the condition in line 2 has 
not been fulfilled then, but is now due to the effect 
of line 5. So we check again:] 
if ngbha\{6} 

remove_V0P (a, 6) 
end if 

[Now we can be sure that the the condition 
(ngbh c\{a,b} = {} or V0P[c] £ Z) is fulfilled for 
c — a,b and we may use the lookup table (cf. 
Eg. 0)j.y 
if {a, b} e E : 

edge +— true 

else: 

edge <— false 
end if 

(edge, V0P[a], VOP [fa]) <- 
cphase.t able [edge , VOP [a] , VOP [b] ] 

remove_VQP (vertex a, vertex b): 

[This reduces VOP[a] to I, avoiding (if possible) to use b 
as swapping partner.] 

[First, we choose a swapping partner c.[ 

if ngbha\{6}/{}: 

c <— any element of ngbha\{£>} 

else: 

c^b 
end if 

d <— decomposition_lookup_table [a] 
[c contains now a decomposition such as Eq. Jfyj 
for v from last factor of d to first factor of d 
ifv = V^iX: 

local_complementation (a) 
else: ( this means that v = \fiZi) 
local_complementation (b) 
end if 
[Now, V0P[a] = I.] 

local_complementation (vertex a) 
[performs the operation specified in Definitional 
n v <— ngbhw 
for i G n v : 

for j £ n v : 
if i < j: 

if e E: 

remove edge (i,j) 

else: 

add edge 
end if 
end if 
end for 

V0P[i] <- V0P[{\V-iZ 
V0P[v] «- V0P[v]Vt3( 
end for 



changes as follows: 



LISTING 1: Pseudo-code for controlled phase gate (AZ) 
acting on vertices a and b (cphase), and for the two 
auxiliary routines remove_VQP and local_complementation. 



\(V,E)) = 



X a I] Zb \ H a \{V,E\{{a,b}\benghha})). 

6^ ngbh a 



(*) 



(12) 



The value of £ is chosen at random (using a pseudo- 
random number generator). To update the simulator 
state, the VOPs are right-multiplied with the under- 
braced operators (*) and the edges incident on a are 
deleted as indicated in the ket. 

A measurement of the Y observable (P — ±F) requires 
a complementation of the edges set according to 

£h E A {{b, c} | b, c S ngbh a} 

and a change in the VOPs as follows: 



Cfc i ► Cb^—iZ ^ for b G ngbh a U {a}, 

where the dagger in parentheses is to be read only for 
measurement result f = 1. 

The most complicated case is the X measurement 
which requires an update of edges and VOPs as follows: 

E h-> E A {{c, d} | c e ngbh b, d £ ngbh a} 
A {{c, d} | c, d <G ngbh b (1 ngbh a} 
A {{b,d} | d£ ngbha\{6}} 



C C Z^ for c = a 
-(t) 



aviv v 



c r z 



for c = 6 (read only for £ = 1) 

ngbh a \ ngbh 6 \ {6} 
(for C = 0) 

ngbh b \ ngbh a \ {a} 
(for C = 1) 



for c £ < 



otherwise 



(13) 



Here, & is a vertex chosen arbitrarily from ngbh a and 
/ 1 -1 \ 



^ ^ i i j 

In all these cases the measurement result is chosen at 
random. Only in case of the measurement of P a = ±X an 
isolated vertex, the result is always £ = (which means 
an actual result of C = for P a — X and £ = 1 for 

Pa = -X.) 
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1 import random 

2 import graphsim 

3 

4 gr = graphsim. GraphRegister (8) 

5 

e gr.hadamard (4) 

7 gr.hadamard (5) 

s gr.hadamard (6) 

o gr.cnot (6, 3) 

io gr.cnot (6, 1) 

n gr.cnot (6, 0) 

12 gr.cnot (5, 3) 

13 gr.cnot (5, 2) 

14 gr.cnot (5, 0) 
is gr . cnot (4 , 3) 

16 gr.cnot (4, 2) 

17 gr.cnot (4, 1) 

18 

19 for i in xrange (7) : 

20 gr.cnot (i, 7) 

21 

22 print gr. measure (7) 

23 

24 gr .print_stabilizer () 

LISTING 2: A simple example in Python 

V. IMPLEMENTATION 

The algorithm described above has been implemented 
in CH — h in object-oriented programming style. We have 
used the GNU Compiler Collection (GCC) 24] under 
Linux, but it should be easyto compile the program 
on other platforms as well 28] . The implementation 
is done as a library to allow for easy integration into 
other projects. We also offer bindings to Python J2l|, so 
that the library can be used by Python programs as well. 
(This was achieved using SWIG |2q.1 

The simulator, called "GraphSim" can be downloaded 
from i|. 

A detailed documentation of the library is supplied 
with it. To demonstrate the usage here at least briefly, 
we give Listing 2 as a simple toy example. It is written 
in Python, and a complete program. 

In the example, we start by loading the GraphSim li- 
brary (Line 2) and then initialize a register of 8 qubits 
(line 4), which are then all in |0) state. We get an object 
called "gr" of class GraphRegister, which represents the 
register of qubits. For all following operations, we use 
the methods of gr to access its functionality. In our ex- 
ample, we simply build up an encoded "0" state in the 
well-known 7-qubit Steane code, which we then measure. 

First, we apply Hadamard and cnot gates onto the 
qubits with number through 6 in order to build up 
the Steane-encoded "0" (Lines 6-17). To check that we 
did so, we measure the encoded qubit, which is done by 
using cnot gates to sum up their parity in the eighth 
qubit ("qubit 7") (Lines 19, 20). Measuring qubit 7 then 
gives "0", as it should (Line 22). 

For further details on using of the GraphSim library 



from a C++ or Python program, please see the docu- 
mentation supplied with the source code [23|. 

With approximately 1400 lines, GraphSim is complex 
enough that one cannot take for granted that it faithfully 
implements the described algorithm without bugs, and 
testing is necessary. Fortunately, this can be done very 
conviniently by comparing with Aaronson and Gottes- 
man's "CHP" simulator. As these two programs use 
quite different algorithms to do the same task, it is very 
unlikely that any bugs, which they might have, produce 
the same false results. Hence, if both programs give the 
same result, they can reasonably be considered both to 
be correct. 

We set up a script to do random gates and measure- 
ments on a set of qubits for millions of iterations. All op- 
erations were performed simultaneously with CHP and 
GraphSim. For measurements whose outcome was cho- 
sen at random by CHP, a facility of GraphSim was used 
that overrides the random choice of measurement out- 
comes and instead uses a supplied value. For measure- 
ments with determined outcome, however, it was checked 
whether both programs output the same result. Also, 
every 1000 steps, the stabilizer tableau of GraphSim's 
state was calculated from its graph representation and 
compared to CHP's tableau. |29j 

After simulation 4 • 10 6 operations on 200 qubits in 18 
hours and 2 • 10 8 operations on 20 qubits in 19.7 hours 
without seeing discrepancies, we are confident that we 
have exhausted all special cases, so that the two programs 
can be assumed to always give the same output. As they 
are based on very different algorithm, this reasonably 
allows to conclude that they both operate correctly. 



VI. PERFORMANCE 

We now show that our simulator yields the promised 
performance, i. e. performs a simulation of M steps in 
time of order O(NdM), where N is the number of qubits 
and d the maximum vertex degree that is encountered 
during the calculation. Let us go through the different 
possible simulation steps in order to assess their respec- 
tive time requirements. 

Single-qubit gates are fastest: they only need one look- 
up in the multiplication table of the local Clifford group 
(which is hard-coded into the simulator), and are hence 
of time complexity 6(1). 

Measurements have a complexity depending on the ba- 
sis in which they have to be carried out. For a Z measure- 
ment, we have to remove the deg a edges of the measured 
vertex a. As d is the maximum vertex degree that is to 
be expected within the studied problem, the complexity 
of a Z measurement is 0{d) < 0(N) (as d < N). 

For a Y and X measurement, we have to do local com- 
plementation, which requires dealing with up to d ^ d 2 ^ 
edges, and hence, the overall complexity of measurements 
is 0(d 2 ). 
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G O CHP 

□ O GraphSim 



size of clusier state ( 1 (Kill stales of ihis size) 

FIG. 2: Comparison of the performance of CHP and Graph- 
Sim. A simulation of entanglement purification was used as 
sample application. The register has 1000 times the size of 
the states to hold an ensemble of 1000 states. 



10000 100000 
width ol register (number of qubits) 



FIG. 3: Benchmark of GraphSim for very large registers. En- 
tanglement purification -specifically: the purification of 10- 
qubit cluster states with the protocol of Ref. 0- was used as 
sample problem. The register was filled up with cluster states 
to make a large ensemble, and two protocol steps were simu- 
lated. The average time per operation was obtained from the 
total run-time 12711. 



For the phase gate, the same holds. Here, we need a 
fixed number (up to 5) of local complementations. Thus, 
measurements and two-qubit gates take 0(d 2 ) time. 

This would be no improvement to Aaronson and 
Gottesman's algorithm, if we had d = O(N). The latter 
is indeed the case if one applies randomly chosen opera- 
tions as we did to demonstrate GraphSim's correctness. 
There, we indeed did not observe any superiority in run- 
time of GraphSim. 

In practice, however, this is quite different. For ex- 
ample, when simulating quantum error correction, one 
can reasonable assume d = (log AT). This is because 
all QEC schemes avoid to do to many operations on one 
and the same qubit in a row, as this would spread errors. 
So, vertex degrees remain small. The same reasoning 



applies to entanglement purification schemes and, more 
generally, to all circuits which are designed to be robust 
against noise. 

The space complexity is dominated by the space 
needed to store the quantum state representation. As 
argued in Section this requires only space of O(Nd), 
where d is the average vertex degree. As explained above, 
we may expect d (as d) to scale sub-linearly with N in 
typical application, in many applications as O (N log N). 
This is what allows us to handly substantially more 
qubits than it is possible with the 0(N 2 ) tableau rep- 
resentation. 

As a first practical test, we used GraphSim to simu- 
late entanglement purification of cluster states with the 
protocol of Ref. 7] . This has been a starting point of a 
detailed analysis of the communication costs of establish- 
ing multipartite entanglement states via noisy channels 
[l3| . Fig. 1^1 demonstrates that GraphSim is indeed suit- 
able for this purpose. Note, that for the right-most data 
points, the register holds 30,000 qubits. 

As we did a Monte Carlo simulation, we had to loop 
the calculation very often and still got an output within 
a few hours. For simulations involving several millions 
of qubits and a large number of runs, we waited about a 
week for the results when using eight processors in par- 
allel. We redid some of these calculations in a more con- 
trolled testing environment as a benchmark for Graph- 
Sim. Fig. shows the results in a log- log plot. 



VII. CONCLUSION 

To summarize, we have used recent results on graph 
states to find a very space-efficient representation of sta- 
bilizer states, and determined, how this representation 
changes under the action of Clifford gates. This can be 
used to simulate stabilizer circuits more efficiently than 
previously possible. The gain is not only in simulation 
speed, but also in the number of manageable qubits. In 
the latter, at least two orders of magnitude are gained. 
We have presented an implementation of our simulation 
algorithm and will soon publish results about entangle- 
ment purification which makes use of our new technique. 
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